{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一.什么是隐马尔可夫模型\n",
    "隐马尔科夫模型（HMM）主要用于时序的概率模型，它假设由一个隐藏的马尔科夫链随机生成不可观测的隐状态序列，再由各个隐状态随机生成一个观测态从而得到可观测序列，如下图：   \n",
    "\n",
    "![avatar](./source/12_HMM定义.png)  \n",
    "\n",
    "所以，HMM可以由初始概率分布，状态转移概率矩阵，观测概率矩阵确定，我们可以对其进行如下定义：   \n",
    "\n",
    "设$Q$是所有可能的状态的集合，$V$是所有可能的观测的集合：   \n",
    "\n",
    "$$\n",
    "Q=\\{q_1,q_2,...,q_N\\},V=\\{v_1,v_2,...,v_M\\}\n",
    "$$  \n",
    "\n",
    "这里，$N$是可能的状态数，$M$是可能的观测数，隐状态序列和观测序列可分别表示为：   \n",
    "\n",
    "$$\n",
    "I=(i_1,i_2,...,i_T),O=(o_1,o_2,...,o_T),i_j\\in Q,o_j\\in V,j=1,2,..,T\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "所以状态转移概率矩阵可以表示为：   \n",
    "\n",
    "$$\n",
    "A=[a_{ij}]_{N\\times N}\n",
    "$$  \n",
    "这里,$a_{ij}=P(i_{t+1}=q_j\\mid i_t=q_i),i=1,2,...,N,j=1,2,...,N$  \n",
    "\n",
    "观测概率矩阵可以定义为：   \n",
    "\n",
    "$$\n",
    "B=[b_j(k)]_{N\\times M}\n",
    "$$  \n",
    "\n",
    "这里，$b_j(k)=P(o_t=v_k\\mid i_t=q_j),k=1,2,..,M,j=1,2,...,N$   \n",
    "\n",
    "初始状态的概率分布，可以定义为：   \n",
    "\n",
    "$$\n",
    "\\pi=(\\pi_i),\\pi_i=P(i_1=q_i),i=1,2,...,N\n",
    "$$   \n",
    "\n",
    "所以,HMM的参数：  \n",
    "\n",
    "$$\n",
    "\\lambda=(A,B,\\pi)\n",
    "$$   \n",
    "\n",
    "HMM的联合概率分布也能很容易写出来：   \n",
    "\n",
    "$$\n",
    "p(o_1,o_2,...,o_T,i_1,i_2,...,i_T)=p(i_1)p(o_1\\mid i_1)p(i_2\\mid i_1)p(o_2\\mid i_2)\\cdots p(i_T\\mid i_{T-1})p(o_T\\mid i_T)\\\\\n",
    "=p(i_1)p(o_1\\mid i_1)\\prod_{t=2}^Tp(i_t\\mid i_{t-1})p(o_t\\mid i_t)\\\\\n",
    "=\\pi_{i_1}b_{i_1}(o_1)\\prod_{t=2}^Ta_{i_{t-1}i_t}b_{i_t}(o_t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二.$P(O\\mid\\lambda)$概率计算\n",
    "\n",
    "我们通常需要处理这样的任务，给定一个隐马尔科夫模型$\\lambda=(A,B,\\pi)$，计算某观测序列$O=(o_1,o_2,...,o_T)$出现的概率，即$P(O\\mid\\lambda)$，但是隐马尔科夫模型能给我们的是联合概率分布$P(O,I\\mid\\lambda)$，这里$I=(i_1,i_2,...,i_T)$，所以我们需要将隐状态积分积掉：   \n",
    "\n",
    "$$\n",
    "P(O\\mid\\lambda)=\\sum_IP(O,I\\mid\\lambda)\\\\\n",
    "=\\sum_{i_1=1}^N\\sum_{i_2=1}^N\\cdots\\sum_{i_T=1}^N\\pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2)\\cdots a_{i_{T-1}i_T}b_{i_T}(o_T)\n",
    "$$   \n",
    "\n",
    "但如果直接这样暴力求解，计算的时间复杂度很高，为$O(TN^T)$，解决方案其实在马尔科夫链的计算中已经有提及了，我们可以采用动态规划的方式，将需要重复计算的结果缓存起来，我们简单分析一下，对于每个求解项，当对于$i_2,i_3,...,i_T$进行求和计算时，其实每项前面的$\\pi_{i_1}b_{i_1}(o_1),i_1=q_1,q_2,...,q_N$都会被重复计算多遍，同理对于$i_3,i_4,...,i_T$进行计算时，包含$i_2,o_2$及其前面的项也会被重复计算...，所以如果在每次计算是都利用前面已经计算好的结果时，这将大大减少计算量，具体说来我们可以这样操作：   \n",
    "\n",
    "#### 前向算法\n",
    "首先对前向概率做定义$\\alpha_t(i)=P(o_1,o_2,...,o_t,i_t=q_i\\mid\\lambda)$,即到时刻$t$，部分观测序列为$o_1,o_2,...,o_t$且隐状态为$q_i$的概率，接下来叙述算法流程：   \n",
    "\n",
    ">（1）计算初始值：   \n",
    "\n",
    "$$\n",
    "\\alpha_1(i)=\\pi_ib_i(o_1),i=1,2,...,N\n",
    "$$   \n",
    "\n",
    "> （2）递推 对$t=1,2,...,T-1$，计算：   \n",
    "\n",
    "$$\n",
    "\\alpha_{t+1}=[\\sum_{j=1}^N\\alpha_t(j)a_{ji}]b_i(o_{t+1}),i=1,2,...,N\n",
    "$$  \n",
    "\n",
    "> （3）最后：   \n",
    "\n",
    "$$\n",
    "P(O\\mid\\lambda)=\\sum_{i=1}^N\\alpha_T(i))\n",
    "$$   \n",
    "\n",
    "可以发现使用动态规划后，复杂度大大降低了，只有$O(N^2T)$，由于条件独立性假设，后面的概率计算其实与前面相互独立，所以我们也可以从反方向使用动态规划   \n",
    "\n",
    "#### 后向算法\n",
    "\n",
    "首先对后向概率做定义$\\beta_t(i)=P(o_{t+1},o_{t+2},...,o_T\\mid i_t=q_i,\\lambda)$，即在时刻$t$状态为$q_i$的条件下，从t+1到T的部分观测序列为$o_{t+1},o_{t+2},...,o_T$的概率，接下来看看求解流程：   \n",
    "\n",
    ">（1）赋初始值：   \n",
    "\n",
    "$$\n",
    "\\beta_T(i)=1,i=1,2,...,N\n",
    "$$   \n",
    "\n",
    ">（2）对$t=T-1,T-2,...,1$：   \n",
    "\n",
    "$$\n",
    "\\beta_t(i)=\\sum_{j=1}^Na_{ij}b_j(o_{t+1})\\beta_{t+1}(j),i=1,2,...,N\n",
    "$$   \n",
    "\n",
    "> （3）最后：   \n",
    "\n",
    "$$\n",
    "P(O\\mid\\lambda)=\\sum_{i=1}^N\\pi_ib_i(o_1)\\beta_1(i)\n",
    "$$    \n",
    "\n",
    "#### 前向+后向\n",
    "\n",
    "其实前向和后向算法是可以同时进行的，它们互相独立计算，然后到中间任意时刻$t,t=1,2,...,T$再“会师”即可，如下图：   \n",
    "\n",
    "![avatar](./source/12_HMM前向后向.png)   \n",
    "\n",
    "所以：   \n",
    "\n",
    "$$\n",
    "P(O\\mid\\lambda)=\\sum_{i=1}^N\\alpha_t(i)\\beta_t(i),t=1,2,...,T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三.代码实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "隐马尔科夫模型实现，封装到ml_models.pgm\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "class HMM(object):\n",
    "    def __init__(self, hidden_status_num=None, visible_status_num=None):\n",
    "        \"\"\"\n",
    "        :param hidden_status_num: 隐状态数\n",
    "        :param visible_status_num: 观测状态数\n",
    "        \"\"\"\n",
    "        self.hidden_status_num = hidden_status_num\n",
    "        self.visible_status_num = visible_status_num\n",
    "        # 定义HMM的参数\n",
    "        self.pi = None  # 初始隐状态概率分布 shape:[hidden_status_num,1]\n",
    "        self.A = None  # 隐状态转移概率矩阵 shape:[hidden_status_num,hidden_status_num]\n",
    "        self.B = None  # 观测状态概率矩阵 shape:[hidden_status_num,visible_status_num]\n",
    "\n",
    "    def predict_joint_visible_prob(self, visible_list=None, forward_type=\"forward\"):\n",
    "        \"\"\"\n",
    "        前向/后向算法计算观测序列出现的概率值\n",
    "        :param visible_list:\n",
    "        :param forward_type:forward前向，backward后向\n",
    "        :return:\n",
    "        \"\"\"\n",
    "        if forward_type == \"forward\":\n",
    "            # 计算初始值\n",
    "            alpha = self.pi * self.B[:, [visible_list[0]]]\n",
    "            # 递推\n",
    "            for step in range(1, len(visible_list)):\n",
    "                alpha = self.A.T.dot(alpha) * self.B[:, [visible_list[step]]]\n",
    "            # 求和\n",
    "            return np.sum(alpha)\n",
    "        else:\n",
    "            # 计算初始值\n",
    "            beta = np.ones_like(self.pi)\n",
    "            # 递推\n",
    "            for step in range(len(visible_list) - 2, -1, -1):\n",
    "                beta = self.A.dot(self.B[:, [visible_list[step + 1]]] * beta)\n",
    "            # 最后一步\n",
    "            return np.sum(self.pi * self.B[:, [visible_list[0]]] * beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "就用李航老师的例10.2做个计算：   \n",
    "\n",
    "（1）开始，假设有3个盒子，选择每个盒子的概率分别为[0.2,0.4,0.4]；   \n",
    "\n",
    "（2）然后，从当前盒子转移到下一个盒子的规则如下,即如果当前盒子是1，那么下一个继续是盒子1的概率为0.2，转移到盒子2的概率是0.2，转移到盒子3的概率是0.3，如果当前是盒子2，转移到盒子1的概率是0.3，继续留在盒子2的概率是0.5...     \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "0.5 & 0.2 & 0.3\\\\ \n",
    "0.3 & 0.5 & 0.2\\\\ \n",
    "0.2 & 0.3 & 0.5\n",
    "\\end{bmatrix}\n",
    "$$   \n",
    "\n",
    "（3）确定盒子后，从每个盒子里面抽出红球、白球的概率矩阵如下，即从盒子1抽出红球的概率为0.5，抽出白球的概率为0.5....   \n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "0.5 & 0.5\\\\ \n",
    "0.4 & 0.6\\\\ \n",
    "0.7 & 0.3\n",
    "\\end{bmatrix}\n",
    "$$   \n",
    "\n",
    "目前观测到的结果为：(红,白,红)，求其概率$P(O\\mid\\lambda)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#造数据\n",
    "pi = np.asarray([[0.2], [0.4], [0.4]])\n",
    "A = np.asarray([[0.5, 0.2, 0.3],\n",
    "                [0.3, 0.5, 0.2],\n",
    "                [0.2, 0.3, 0.5]])\n",
    "B = np.asarray([[0.5, 0.5],\n",
    "                [0.4, 0.6],\n",
    "                [0.7, 0.3]])\n",
    "\n",
    "hmm = HMM()\n",
    "hmm.pi = pi\n",
    "hmm.A = A\n",
    "hmm.B = B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.130218"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#查看结果:前向\n",
    "hmm.predict_joint_visible_prob([0, 1, 0],forward_type=\"forward\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.130218"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#查看结果:后向\n",
    "hmm.predict_joint_visible_prob([0, 1, 0],forward_type=\"backward\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
